Path integral Monte Carlo simulation of the second layer of 4 He adsorbed on graphite 
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We have developed a path integral Monte Carlo method for simulating helium films and apply it 
to the second layer of helium adsorbed on graphite. We use helium-helium and helium-graphite 
interactions that are found from potentials which realistically describe the interatomic interactions. 
The Monte Carlo sampling is over both particle positions and permutations of particle labels. From 
the particle configurations and static structure factor calculations, we find that this layer possesses, 
in order of increasing density, a superfluid liquid phase, a y/7 x \/7 commensurate solid phase that 
is registered with respect to the first layer, and an incommensurate solid phases. By applying the 
Maxwell construction to the dependence of the low-temperature total energy on the coverage, we 
are able to identify coexistence regions between the phases. From these, we deduce an effectively 
zero-temperature phase diagram. Our phase boundaries are in agreement with heat capacity and 
torsional oscillator measurements, and demonstrate that the experimentally observed disruption 
of the superfluid phase is caused by the growth of the commensurate phase. We further observe 
that the superfluid phase has a transition temperature consistent with the two-dimensional value. 
Promotion to the third layer occurs for densities above 0.212 atom/A 2 , in good agreement with 
experiment. Finally, we calculate the specific heat for each phase and obtain peaks at temperatures 
in general agreement with experiment. 



PACS numbers 67.70.+n, 67.40 Kh 
I. INTRODUCTION 

Helium adsorbed on graphite provides an excellent re- 
alization of a number of nearly two-dimensional (2D) 
phenomena. The helium film grows in a succession of 
distinct, atomically thin layers as the density of the 
adsorbate increases, and as many as seven such, layers 
may be observed on a well-prepared substrateE Con- 
sequently, it is possible to investigate the evolution of 
each layer's phase diagram. A number of experimen- 
tal methods have been usecUfpjj-ithis purpose, includin 
specific heat measurements, Bi§EJ jMptron scattering,ET 
torsional oscillator measurementsJijH and third sound. 
The phase diagrams of the layers nearest the substrate 
are rich. Evidence has been found for self-bound fluid 
phases that are superfluid at low temperatures, a vari- 
ety of registered solid structures, and incommensurate 
solid phases. These phases and the coexistence regions 
that separate them are governed by a delicate balance 
of quantum effects, such as large zero-point motion and 
particle permutations, with adatom and substrate inter- 
actions. 

Much of the early experimental work on the helium- 
graphite system concentrated on the first adaoiJaed layer. 
Several reviews of this work are available Bert 2 ] On the 
other hand, until recently, relatively little information 
was available on the phases of the second and higher lay- 
ers. This situation has changed dramatically over the lasi. 
several years. Extensive heat capacity measurementsu'u 
of the first six layers have been performed, and super- 
fluidity in the higtoj layers has been detected by both 
torsional oscillatorEjIiil and third sound measurements. El 
Taken together, these experiments indicate that the sec- 



ond layer has a unique phase diagram, with superfluid, 
commensurate solid, and incommensurate solid phases. 
No other layer exhibits all three phases. 

Motivated by these experiments, we have undertaken 
a path integral Monte Carlo (P1MC) simulation of 
the second adsorbed layer. We identify a liquid (L) 
phase with an equilibrium density of 0.1750 atom/ A 2 , 
a \fl x spj commensurate triangular solid (C) at 0.1996 
atom/A 2 , and an incommensurate triangular solid (IC) 
phase for densities above 0.2083 atom/A 2 . All cov- 
erage values are for the total adsorbed film. Using 
the Maxwell construction, we determine coexistence re- 
gions between these phases, namely the gas-liquid (G- 
L), liquid-commensurate solid (L-C), and commensurate- 
incommensurate solid (C-IC) phases, at effectively zero 
temperature. Our calculated phase diagram confirms 
the idea that the superfluid phase is mlcirapted by the 
formation of the commensurate solid Jij'LIru We further 
show that the liquid phase behaves like a typical two- 
dimensional superfluid. We also calculate the specific 
heat for each phase and find peaks in general agreement 
with the experimental values. Finally, we observe promo- 
tion to the third layer at a coverage in good agreement 
with experiment. A preliminary repeal of some of our 
findings has been published elsewhere.Eij The present pa- 
per expands and extends this work. 

This paper is arranged in the following manner. Sec- 
provides an overview of what is kn own about 
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the second layer from experiments, 
view previous simulations of helium films and the re- 
lated simulation of two-dimensional helium. We point 
out that none of these simulations, while interesting in 
their own right, exhibit all the phenomena observed in 
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the second-layer phase diagram. Section [Ly presents the 
details of our simulation method, which includes particle 
permutations and realistic particle-particle and particle- 
substrate interactions. The results of our calculations are 
presented in Sec. H. We demonstrate the existence of 
each phase, explain the construction of the second-layer 
phase diagram, and present calculations of properties for 
each phase. 

A. Experimental Overview 

Specific heat measurements have formed the basis for 
constructing the first helium layer's phase diagram, but 
until recently, relatively little work was doncs-on the sec- 
ond layer, with a couple of exceptions. BretzE-3 examined 
this layer under compression of the third and obtained 
evidence for the melting of the incommensurate second 
layer solid. The low density_tange of this layer was ex- 
plored by Polanco and Bretz.E.3 They determined that the 
compression of the first layer by the growth of the sec- 
ond must be taken into account in order to determine the 
phases at low second-layer densities. They interpreted 
their results to indicate that the second layer has gas- 
liquid coexistence at low coverages. 

The heat capacity measurements of Greywall and 
Busch provide the most extensive investigation of the 
second-layer phase diagram. They find evidence for four 
phases: gas, liquid, commensurate solid, and incommen- 
surate solid. These phases are identified in the following 
manner. At low densities, a low, rounded peak occurs 
in the heat capacity. This—has previously been associ- 
ated with the liquid phaseO At low temperatures, the 
heat capacity depends linearly on density roughly be- 
tween 0.13 and 0.16 .atom/A 2 , which is a requirement 
for coexisting phases,ta Thus this region can be identi- 
fied as a gas-liquid coexistence region, with the gas phase 
having negligible density at the lowest temperatures. Ev- 
idence for liquid-commensurate solid coexistence can be 
found between 0.187 and 0.197 atom/A 2 . In this region, 
in addition to the low peak associated with the liquid 
phase, another, larger peak at a higher temperature can 
be observed. The location of the larger peak is inde- 
pendent of coverage, suggesting that it may be associ- 
ated with the melting of a commensurate solid phase. 
Greywall suggested^ that this phase corresponds to the 
commensurate structure proposed earlier for 
3 He on graphite.EZrtL3 A third coexistence region occurs 
between 0.2030 and 0.2080 atom/A 2 , where the commen- 
surate melting peak is accompanied by another, lower 
temperature peak. This second peak is associated with 
the melting of an incommensurate solid phase. For cov- 
erages from 0.2080 to the beginning of third layer promo- 
tion at 0.212 atom/A 2 , the incommensurate melting peak 
is the sole feature in the specific heat measurements. Un- 
like the peak associated with the commensurate phase, 
the incommensurate melting peak is temperature depen- 



dent, occurring at about 1 K at the lowest incommen- 
surate densities, but increasing to about 1.5 K at the 
density where third layer promotion begins. 

The principal limitation on using the heat capacity 
measurements to determine the phase diagram is that 
they can only identify phases indirectly so additional 
confirmation is desirable. Direct evidence for the income 
mensurate solid phase comes from neutron scattering ,IJlI 
but no similar evidence exists for the commensurate 
phase. Apparently, the incommensurate phase can be 
resolved in these experiments only after some additional 
compression by the third layer. Consequently, there is no 
scattering evidence for the commensurate solid, which is 
replaced by the incommensurate solid before promotion 
to the third layer begins. 

Further insight into the second-layer phase diagram 
comes from the torsional oscillator measurements of 
Crowell and ReppyEM They detected superfluidity at 
intermediate densities, which incidentally provided direct 
evidence that the second layer has a liquid phase. Ques- 
tions remain about the liquid phase, however, since the 
apparent onset density is somewhat higher than would be 
expected from either the heat capacity measurements or 
the liquid equilibrium density of purely two-dimensional 
heliumJla The range of superfluid coverage also provides 
additional, although indirect, evidence that a solid phase 
begins to form above 0.187 atom/A 2 . Above this density, 
the superfluid signal vanishes and does not reappear un- 
til the third layer. This disappearance coincides almost 
exactly with the growth of the supposed commensurate 
solid phase. Apparently, the growing solid phase disrupts 
the connectivity required to detect superfluidity. 

B. Previous Simulations 

The results of Monte Carlo calculations are often 
used to help interpret the experimental results discussed 
above. The simplest way to treat a helium layer is as 
a purely two-dimensional system, for which there are 
both zero temperature.-and finite temperature calcula- 
tions. Whitlock et al.Ej used Green's function Monte 
Carlo to calculate the equilibrium liquid coverage at zero 
temperature, obtaining 0.04356 atom/A 2 . They also de- 
termined that 2D helium would solidify, and that the liq- 
uid and solid phases coexisted between 0.0678 and-0.0721 
atom/ A 2 . More recently, Gordillo and CeperleyES have 
investigated the 2D phase diagram at finite temperatures 
with path integral Monte Carlo. Their low temperature 
results are consistent with the zero temperature calcula- 
tions. They also determined spinodal lines and found a 
finite density gas phase at temperatures above 0.75 K. 
The direct comparison of these results with the second 
helium layer is limited, since the 2D calculations do not 
include any substrate features and do not allow the film 
to spread perpendicularly. As a result, no commensurate 
solid phase or layer promotion can occur. 
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Simulations of helium films using realistic models for 
the graphite substrate have also been performed. Abra- 
ham and BroughtonEil used path integral Monte Carlo to 
investigate the first layer of 3 Hc on graphite. They were 
able to identify fluid, commensurate solid, domain wall 
liquid and solid, and incommensurate solid phases. No- 
tably, they determined that particle permutations were 
unimportant for the first layer for the coverages they in- 
vestigated, so there was no possibility for superfluidity in 
the simulation. Experimentally, the phase diagrams for 
3 He and 4 He at the densities they simulated arc nearly 
identical, so it reasonable to conclude that their Simula^ 
tion results also apply to 4 He. This work was extendcdta 
to a simulation of the second adsorbed layer of 3 He at 
the v7 x v7 commensurate density. Particle permuta- 
tions were again neglected. It was established that the 
second-layer commensurate phase was stable for-temper- 



atures below 1 K. Very recently, Whitlock et al.crl inves- 
tigated the ground state properties of the first helium 
layer using a laterally averaged potential for the helium- 
graphite interaction. They determined the equilibrium 
liquid coverage and the onset coverage for solidification 
in the first layer, and determined the coexistence region 
between these two phases. They also estimated comple- 
tion densities for the first and second layers, obtaining 
agreement with the experimental results. They did not 
take the corrugations of the graphite substrate into ac- 
count and so did not observe the x -\/3 commensurate 
solid phase that occurs in the first layer. 

Complementary to the ca^ulations discussed above is 
the work by Clements et al£j~E3 using the hypernetted- 
chain Euler-Lagrange theory. For 2D -helium, this 
method reproduces the Monte Carlo resultslij for the liq- 
uid phase and provides a direct calculation of the chem- 
ical potential, third sound, and spinodal points. When 
applied to layered systems, the theory gives liquid cover- 
age ranges and layering transitions but is not capable of 
investigating solid phases. For this reason, these calcula- 
tions are restricted to the third and higher helium layers, 
and assume that the first two layers form an inert, fea- 
tureless solid. Also complementary are the path integral 
Monte Carlo calculations of Wagner and Ceperleyc^O for 
4 He and hydrogen films on crystalline hydrogen. In their 
helium film simulation, superfluidity and layer-by-layer 
growth occurred, but the film did not solidify. 

As we discussed in Sec. I A, the second layer of 4 He on 



graphite is unusual in that it is known experimentally to 
have both a superfluid liquid and two solid phases, one 
commensurate and the other incommensurate with the 
first layer. The simulations discussed above are interest- 
ing in their own right, but none have exhibited the three 
phases seen in the second layer. In order for a simulation 
to produce these phases, it must possess three features. 
First, the presence of superfluidity means that particle 
permutations must be included in the simulation. This 
is because superfluidity results from permutation cycles 
of infinite lengthrS It is also expected that the bound- 
aries of the phases will be effected by permutations. Sec- 



ond, the commensurate second-layer solid is found to be 
registered with respect to the first layer, so the effect of 
first-layer atoms must be taken into account. Third, the 
attraction of the substrate and first layer on the second 
must be implemented correctly so that the commensu- 
rate phase is replaced by the incommensurate phase be- 
fore promotion to the third layer begins. In the following 
section we outline our simulation method, which contains 
the necessary features to exhibit these three phases. 



II. SIMULATION METHODS AND DETAILS 

Path integral Monte Carlo is a powerful tool for simu- 
lating quantum systems at finite temperatures. By incor- 
porating sampling of particle configurations and particle 
permutations, both normal and superfluid helium can be 
simulated. E£l If a substrate is added to the simulation, a 
quantum film will result. The purpose of this section is to 
describe the modifications that are necessary to add the 
effects of the substrate into the simulation. The result 
will be a simulation method that is capable of exhibit- 
ing superfluid, commensurate solid, and incommensurate 
solid phases, as well as layer promotion. 

Central to our PIMC method is the approximation 
used for the high temperature density matrix. It is es- 
sential that the starting temperature be made as low as 
possible so that permutations will be accepted. As we 
will discuss in this section, the graphite substrate com- 
plicates a straight-forward extension of the starting ap- 
proximation used in bulk simulations. For this reason we 
will not include sampling of the first-layer atom configu- 
rations in the calculation and will concentrate instead on 
the second layer. 

It is essential to include the effect of the first layer 
on the second, however. We approximate this effect 
by placing first-layer atoms on the sites of a triangu- 
lar lattice at a fixed height above the substrate. This 
allows us to treat the helium-graphite correlations in 
a much simpler manner, since the atoms on the sec- 
ond layer are not effected by the corrugations of the 
graphite substrate. By not sampling first-layer config- 
urations, we are also able to increase the number of 
second-layer atoms in the simulation. In turn, this allows 
us to scan second-layer coverages in a sufficiently fine 
grid to observe coexistence regions. Having a fine grid 
is particularly important for high second-layer densities, 
since the liquid-commensurate solid and commensurate- 
incommensurate solid coexistence regions exist over rel- 
atively narrow ranges. 

The trade-off for using this approach is that we ignore 
zero-point motion in the first layer. This will cause the 
second layer to form closer to the first layer and have a 
narrower density profiled Ignoring the response of the 
first layer to the second is also known to lead to a lower- 
ing of the [energy of a layer of helium adsorbed onto solid 
hydrogen.cJ However, experimental results indicate that 
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neglecting zero-point motion in the first layer of helium 
on graphite atoms is a reasonable approximation. First, 
the Debye temperature of the solid first layer is greater 
than 50 fCy-and it may be treated as a 2D Debye solid 
up to 3 In our simulation, the temperature is as low 
as 200 mK, and never exceeds 2.2 K, so the first layer is 
relatively stiff. Second, although the first layer is known 
to be compressed by the growing second-layer, this is 
most important at low second-layer densities, just after 
second-layer promotion, begins .E£l The coverages studied 
by Polanco and BretzEj are below the range of our sim- 
ulation. As we shall see, our approach is sufficient to 
reproduce many of the observed features of the second 
layer. 

A. Path integral representation of the partition 
function 

We wish to study the problem of a quantum N-particle 
system in the presence of a substrate. The Hamiltonian 
for this system may be written as 

N N 

i—l i<j 

N 

+ ^2v S ub(ri), (1) 

where V2B is the spherically symmetric two-body poten- 
tial between particles, and u stl & is the external field pro- 
duced by the substrate. The two-body potential for he^ 
hum is accurately represented by the Aziz potential.E£l 
Previous path integral simulations using this potential 
have proven quite capable pfjreproducing numerous prop- 
erties of liquid heliumE3~E3 The potential between he- 
lium r^pd graphite has been investigated by Carlos and 
Cole.Ej Using helium-scattering data, they evaluated sev- 
eral forms for the helium-graphite potential. In order to 
write this potential in a pair form, anisotropic terms that 
effectively enhance corrugation must be included. Of the 
potentials examined, an anisotropic 6-12 Lennard-Jones 
potential was found to be preferable, although the form 
was not uniquely determined. For helium atoms more 
than 4 A above the substrate, corrugations are negligible, 
and the anisotropic potential can be replaced by a later- 
ally averaged potential that depends only on the height 
of the atom above the substrate. 

The statistical mechanics of quantum systems are gov- 
erned by the density matrix and the partition function. 
For a system of N bosons at an inverse temperature (3, 
the density matrix is given by 

p(R, R ; /?) = -L £ < Rle-^IPR' >, (2) 

p 

where R and R' are two configurations of N bosons. The 
sum over P is over all permutations of particle labels, 



and PR' is one such permutation. Permutations lead di- 
rectly to the off-diagonal long-range order that produces 
superfluidity. The partition function, Z, is found by in- 
tegrating the diagonal elements of the density matrix, 

Z = /vTE [p(K,PR,[3)d 3 R- (3) 
■ P J 

Evaluating the partition function for interacting sys- 
tems at very low temperatures is complicated by the fact 
that the kinetic and potential terms in the exponent of 
the density matrix cannot be separated, so the form of 
the density matrix is not known in, for instance, the con- 
figuration space representation. We can avoid this prob- 
lem by inserting M — 1 intermediate configurations into 
Eq. (H) to obtain the path integral formulation of the 
partition function, 

Z = 7vtE/-/ d 3 Ri...d 3 R M _ 1 d 3 R 

xp(R, R 1 ;t)/;(R 1 ,R 2 ;t) . . . p(R M -i, PR; t), (4) 

where r = (3/M. The problem of evaluating the partition 
function at a low temperature, /3 _1 , has been replaced by 
the problem of multiple integrations of density matrices 
at a higher temperature, r _1 . The advantage of this is 
that the high temperature density matrices may be ap- 
proximated. In practice, the integrals appearing in Eq. 
([|) cannot be directly evaluated for systems of strongly 
interacting particles. Monte Carlo sampling may be used 
instead to generate configurations and calculate observ- 
ables. 

Equation (Q) lends itself to an interesting visualization. 
The N quantum particles can be thought of as N inter- 
acting classical ring polymers, each with M beads. Sam- 
pling the partition function then corresponds to sampling 
the possible configurations of these polymers. Further- 
more, particle permutations may be introduced into the 
Monte Carlo method by splicing together two or more 
polymer chains. This is known as the polymer isomor- 
phism. 

B. Approximating the density matrix 

In order to use Monte Carlo sampling on the partition 
function, we must first provide a starting approximation 
for the high temperature density matrices that appear in 
the integrand of Eq. (|4|). The simplest starting approxi- 
mation is to use a very large M, which allows us to sepa- 
rate the density matrix into kinetic and potential energy 
terms. This is the semiclassical approximation and is 
exact in the limit M — > oo, according to the Trotter the- 
orem. For superfluid helium systems it is necessary to go 
beyond the semiclassical approximation so that the start- 
ing temperature may be made as low as possible. This 
makes sampling the permutations feasible and speeds the 
equilibration of the ring polymers by avoiding excessively 
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long chains. The high-temperature density matrix we in- 
troduce below can be used with starting temperatures 
as low as 40 K. We thus only have to use, for instance, 
M = 40 to simulate a system at 1 K. 

We approximate the high temperature density matrix 
as a product of the exact free particle solution, an effec- 
tive two-body interaction found from the exact solution 
for two interacting helium atoms, and an effective exter- 
nal interaction found from the exact solution for a single 
atom in a graphite potential: 

N 

p(R,R';r) « n^ ree ( r *> r - T ) 

i=i 

N 

i=l 
N 

x n^ e ( r «' r 'y ;T )' ( 5 ) 

i<j 

where r.^ = — rj. The terms pf ree , p^ e , and p Gr 
will be discussed below. This approximation assumes 
that three-body contributions are negligible and that the 
helium-helium and helium-graphite interactions can be 
decoupled. The former has been shown to be valid for 
bulk helium systems with starting temperatures as low 
as 40 K. 

The term p{ ree is the density matrix for a free particle 
of mass m, given by 



free 



(r,r';r) = A t - 3 cxp[- 7 r(r-r') 2 A?]. 



(6) 



where A t = \J2ittTi 2 /m is the mean thermal wavelength 

for the temperature f/r. 

The helium-helium term, pf^ e , is the interacting part of 
the solution to the density matrix for two helium atoms. 
This can be found by separating the density matrix into 
center-of-mass and relative components. The density ma- 
trix for the relative coordinates is a solution to 



— (r ij) r , ij;T )= [ (h 2 /m)W 2 

-V^in^in^r'^r). (7) 

This equation is equivalent to that satisfied by the time 
evolution propagator in imaginary time. We solve this, 
equation using the methods discussed by Ceperley.E£l 
Briefly, the density matrix can be expanded in a se- 
ries of partial waves and the expansion coefficients are 
found by using the matrix-squaring method. The re- 
sulting solution is used to define the effective helium- 
helium interaction, U (r^ , r'ij ; r) = — \n(p He ) where 
p He = p He /p' ree . This is a six-dimensional function, 
but the spherical symmetry of the density matrix allows 
us to approximate it as a series of one-dimensional func- 
tions. This greatly reduces the memory requirements and 
increases the speed at which the density matrix can be 
evaluated for a particular configuration. 



The density matrix for a single helium atom above a 
graphite substrate is a solution to 



Gr 



-(r,r';r) = [(h 2 /2m)V 2 - V Gr (r)]p Gr (r,r< ;r), (8) 



where V Gr (r) is the full graphite potential. The helium- 
graphite term, p Gr , is the interacting part of the solution 
to this equation. Near the substrate, the potential V 
is anisotropic. A straight forward solution to Eq. (||) 
is to solve it at grid points within a graphite unit cell 
using, for instance, a three-dimensional implicit method 
with periodic boundary conditions at the edges of the 
cell. The resulting six-dimensional function can be ap- 
proximated as a series, expanding around the diagonal 
elements, but this still gives a series of three-dimensional 
functions. This greatly complicates Monte Carlo simu- 
lations of the first-layer atoms using Eq. (||), since stor- 
age requirements become large and evaluating the den- 
sity matrix by interpolating from three-dimensional ta- 
bles becomes excessively burdensome. Thus, simulating 
the first adsorbed layer using a high-temperature density 
matrix is a much more complicated problem than simu- 
lating bulk helium. One could always avoid these prob- 
lems by simply starting at a high enough temperature so 
that the semiclassical approximationcil can be used for 
atoms near the substrate, but then getting permutations 
accepted becomes exceedingly unlikely. 

The problem becomes much simpler further above the 
substrate, where corrugations may be ignored. The 
helium-graphite potential can be found by laterally av- 
eraging over the surface, eliminating the x-y plane peri- 
odicity that complicates the solution near the substrate. 
The helium atom experiences only a z-dependent poten- 
tial, so Eq. (||) can be solved by separating pf r (r,r',r) 
into x, y and z components. The x and y components 
are one-dimensional, free-particle density matrices. The 
solution for p(x,x';t), for instance, is 



pf ree (x, x'; r) = Aj" 1 exp[-^(x - x') 2 /X 2 



(9) 



A similar solution exists for p(y, y'; r). The z-dependence 
is found by solving the parabolic partial differential equa- 
tion 



dp 
#7 



(z, z'; t) = [(h 2 /2m)d 2 /dz 2 - V Gr (z)]p(z, z'; r) 



(10) 

where V Gr (z) is the laterally averaged potential.!^ This 
can be .solved by matrix squaring, or by an implicit 
method.tll The initial condition is that the density ma- 
trix is a delta function at r = 0. We define the ef- 
fective interaction for the helium-graphite density ma- 
trix, U Gr (z,z';T) = -\n[p(z,z';T)/pf ree (z,z'- 1 T)}. This 
is still a function of two variables. In order to make 
evaluating the density matrix efficient during the Monte 
Carlo runs, we expand U as a series of one dimensional 
functions. We rewrite U Gr (z,z') = U(z,Az), where 
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z = (z + z')/2 and Az — \z — z'\. The matrix is dom- 
inated by the diagonal elements, so we expand it as a 
series about (Az) 2 : 



U<* r (z,z',r) = 



U° r (z,Z,T) + U Gr ( Z ',z',T) 

2 

J2u m (z)(^) 2m - 



(11) 



2.0 



1.0 



0.0 



-i.o -\, 



+ Exact 

Series 

Endpoinl 



The average over the two diagonal parts of the solu- 
tion in the first term is called the endpoint approxima- 
tion. The functions U m (z) are found by \ 2 fitting Eq. 
( pi] ) to the exact solution. One simply terminates the 
series when the approximation is sufficiently close to the 
exact solution. Results for the diagonal solution and the 
first two expansion terms are shown in Fig. The off- 
diagonal terms become negligible for z > 4A. The di- 
agonal solution can be compared with the semiclassical 
approximation. Figure 2 compares the exact solution for 
off-diagonal elements to the expansion, Eq. (O), and the 
endpoint approximation, 1/2[tV(z) + tV(z'J[7 
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FIG. 2. The exact solution U Gr {z SL z', t) for z' = 2.821 
compared with the expansion, Eq. (|ll|), and the endpoint 
approximation using tV. 
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FIG. 1. The diagonal and lowest-order off-diagonal terms 
of the expansion of U Gr , Eq. ([n]). The semiclassical approx- 
imation is also shown. The laterally averaged potential was 
used and r = 0.025A'- 1 . 



C. Sampling the density matrix 

With the first layer frozen, the density matrix, Eq. (||), 
for the active second layer atoms can be written in the 
form p = exp(— S), where 



S(R,R';t) = (3iV oct /2)ln(A 2 ) 

N act Naet 

2 



2 , , tt(R-R') 2 



A? 



where = |r 7 ; 



i=l J=l 

N act 

+ ^ ( 7 Gr (^,^;r), (12) 

i=l 

j | . The number of active and frozen 
helium atoms is given by N act and Nj r , respectively. In 
the polymer isomorphism, S is the action for a system 
of interacting polymers. In sampling the paths, we are 
effectively choosing between two different polymer con- 
figurations. The one with the lower action is the more 
favorable configuration, and is more likely to be chosen 
in a Metropolis-style acceptance test. 

As in standard Monte Carlo simulations, the inter- 
action U He is cut off at some maximum distance r c < 
min(L.j;, L y ), where L x and L y are the dimensions of the 
simulation cell. The long-range correction to the inter- 
action felt by each particle is, in cylindrical coordinates 
0, z), 



U 



He 
LR 



2tt 



n(z')dz' / pdpU He (r,r;T) 



(13) 



where r = yj p 2 + (z — z') 7 



(z — z') 2 , and only 



diagonal elements need to be considered. The integral 
of n(z') gives the density of the system. We make the 
approximation that the layer thicknesses can be treated 
as delta functions. This is exact for the frozen first layer. 
Then n(z') = n fr S(z' - z fr ) + n act S(z' - z act ) and 
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where n f r and n act are the densities of the first (frozen) 
and second (active) layers. The factor of one-half before 
the contribution from the active layer is needed to avoid 
double counting. A similar long-range correction is added 
to dU He /dr in the energy calculation. 

As we have emphasized, particle permutations must 
be included in simulations of superfluid helium. These 
permutations correspond to splicing together two or more 
of the polymer rings. This splicing can be accomplished 
by proposing cyclic permutations involving one to four 
particle labels on inverse-temperature slice i + n relative 
to slice i, where n — 2 l and I is the overall level of the 
move. The paths followed by the permuted particles on 
the intermediate slices i + 1 to i + n— 1 that produce the 
permutation are then filled in by successively bisecting 
the interval i to i + n. This is known as multilevel Monte 
Carlo sampling, an extension of the standard Metropolis 
method. The interested reader is referred to a recent 
review article on the subject SB 

In our Monte Carlo runs for helium films, we take / = 3, 
since this gives the best balance between accepting sin- 
gle particle and multiple particle moves. Increasing I in- 
creases the number of permutations that can be accepted 
but decreases the overall acceptance rate, while decreas- 
ing I has the opposite effect. The overall acceptance rate 
for the moves varies between 8% and 15%, depending 
on the density. Tests using I = 2 at selected densities 
showed that the I = 3 results had converged. The accep- 
tance rate of multiparticle permutations is small, between 
0.2% and 0.3% in the liquid phase. We have found that 
similar small acceptance rates are sufficient to obtain the 
superfluid density in bulk simulations. 



D. Calculating observables 

The expectation value of an observable, A, can be 
found from the trace, < A >= Z~ x TrAp. We use PIMC 
to calculate expectation values for the total energy, the 
superfluid density, and the static structure factor. Be- 
low we give formulas for each of these calculations for a 
helium film on a substrate. 

The total energy is given by the expectation value 
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AR is the change in the particle positions between two 
consecutive inverse-temperature slices. The terms U^ al 
and V^J tal are shorthand for the sums over the interaction 
terms in Eq. ( |l2| ) . Notice that the zero of the total energy 
occurs at zero second-layer coverage, where there are no 
active atoms. 

The superfluid density can be calculated using the 
winding number, W, for simulations that have periodic 
boundary conditions. Nonzero winding numbers occur 
when particles, through a series of permutations, are per- 
muted with periodic images of themselves. The winding 



number is directly related to p s , the superfluid densityS 
For a system with periodic boundary conditions in the 
x-y plane, the superfluid density is given by 



£a _ m < (W ■ L) 2 > 

P ~ 2f3fi 2 N act ' 
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where the elements L x and L y are the dimensions of the 
simulation cell. 

Finally, structural information can be obtained with 
the static structure factor, 
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We take z to be perpendicular to the plane of the sub- 
strate, so k = (k x ,k y ). p(k) — 53i=i* exp(ik • r^) is the 
Fourier transform of the density. 



E. Testing the method 

As can be seen from the previous discussion, simulat- 
ing helium systems below the superfluid transition is an 
extremely complicated task, and it is important to verify 
all parts of the method. We have verified our calculations 
for the solution to Eq. (Q) by comparing our results to 
published results for the Lennard-JonesE3 potential and 
to the Aziz potential. The solution to Eq. ( [To| ) for the 
helium-graphite density matrix was checked by compar- 
ing the results obtained from the matrix squaring and 
implicit solution methods. We have verified that the full 
Monte Carlo method outlined above works for bulk he- 
lium systems by reproducing reported va luea for the en- 
ergy, specific heat, and superfluid densityOtij We believe 
these tests sufficiently prove that our simulation method 
works and can be extended to helium films. 



F. Choosing simulation cells 

We perform calculations with a variety of simulation 
cells that are appropriate for examining different regions 
of the second-layer phase diagram. The first consider- 
ation is to choose a simulation cell that will match the 
periodicity of the first-layer triangular solid. This can be 
done by using a rectangular unit cell with a two-point 
basis, with unit vectors ai = ax and a 2 = \/3ay, where 
a = 3.015A Two first-layer helium atoms are located in 
each unit cell at hi = and b 2 = a 1 /2 + a 2 /2. This gives 
a coverage of A. 1270 atom/A 2 , the fully compressed first- 
layer densityfl In examining the second layer, our first 
goal is to scan the layer at intermediate and higher den- 
sities by varying the number of particles and to calculate 
the total energy at each density. For these calculations 
we use simulation cells with dimensions (5ai,3a2) and 
(8ai, 5a2), hereafter referred to as the 5x3 cell and the 
8x5 cell, respectively. The number of active particles in 
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calculations using the 5x3 cell ranged from 8 to 21, cor- 
responding to densities 0.1605 to 0.2159 atom/A 2 . Cal- 
culations with the 8x5 cell had 24 to 52 active particles, 
corresponding to densities between 0.1651 and 0.2096 
atom/ A 2 . These two simulation cells are nearly square, 
which is useful for calculating winding numbers. As will 
be discussed in Sec. HI, the energy calculations for the 
5x3 cell are used to verify that finite-size effects are not 
important in the 8x5 cell. Our conclusions about the 
coverage ranges of various phases are drawn from results 
using the 8x5 cell. 

At high second-layer densities, commensurate and in- 
commensurate triangular solid phases occur. In order to 
further investigate these phases, we use simulation cells 
that can contain an integer number of unit cells of both 
the first- and second-layer solids. That is, the simulation 
cells have the periodicity of both the first- and second- 
layer solids. It is also important to note that the solid 
phases will tend to align with the x and y axes of the 
simulation cell. For the incommensurate solid we use a 
cell with dimensions (5a!,5a 2 ), hereafter referred to as 
the 5x5 cell. This cell can accommodate 32 second 
layer atoms in an equilateral triangular lattice. A dia- 
gram of a second-layer incommensurate solid in the 5x5 
cell is shown in Fig. ||. The second-layer solid is in- 
commensurate with respect to the first since no supercell 
with dimensions less than the minimum dimension of the 
simulation box can be drawn in which both first- and 
second-layer atoms are periodically repeated. 




The simulation of the \/7 x y/7 triangular commensu- 
rate solid presents an additional problem since this struc- 
ture is rotated with respect to the first layer. This tri- 
angular solid can be regarded as having a rectangular 
unit cell with a fourteen point basis. The unit vectors 
for this solid are si = 2ai + b 2 and s 2 = — 2ai + a 2 + b 2 . 
Note that |s 2 | = \/3|si| and |sj| = a/7 | a i | , i — 1,2. 
We use simulation cells with dimensions (2si,2s 2 ) and 
(3si, 2s 2 ) to identify the solid configuration and calculate 
the static structure factor. The commensurate density 
0.1996 atom/ A 2 corresponds to 32 and 48 active parti- 
cles, respectively, for these two cells. 




3i 



FIG. 3. Diagram of the 5x5 simulation cell. The shaded 
circles denote positions of the first layer atoms. The 32 open 
circles denote possible positions of atoms in the second-layer 
incommensurate triangular solid. The arrows indicate the 
unit vectors for the solid described in the text. The lines 
emphasize the triangular structure of the solid. 



FIG. 4. Diagram of a simulation cell used for the v7 x v7 
solid. The dimensions are (3si,2s2). The shaded circles de- 
note positions of the first layer triangular solid. The open 
circles denote possible positions of the second layer regis- 
tered solid. The arrows indicate the unit vectors for the solid 
described in the text. The lines emphasize the triangular 
structure of the solid. The heavily shaded lines indicate the 
VT x V7 supercells. 

A diagram of the (3si,2s 2 ) simulation cell with the 
second layer atoms in v7 x y7 registry is shown in Fig. 
[|. The large, rotated rectangle gives the bounds of the 
simulation cell. First layer atom positions outside this 
rectangle are periodic images of interior atoms. Note that 
the location of the origin is arbitrary. It is not necessary, 
for instance, to place it at a high symmetry point of the 
first-layer lattice, such as over a first-layer atom or at a 
potential minimum. The essential requirements for the 
existence of the partially registered solid are that once 
the origin is chosen, all of the two-dimensional space can 
be divided up into periodically repeated superlattice unit 
cells (supercells), and that the relationships of the first- 
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and second-layer atoms to each other and to the super- 
cell are the same in every supercell. We have chosen the 
origin so that the second-layer atoms can be used to di- 
vide up the rectangle into supercclls. These (primitive) 
supercells are the equilateral parallelograms formed by 
the heavily shaded lines in the interior of Fig. f|. They 
can be seen to exactly fill the rectangle. Second layer 
atoms are located at the four corners, on each of the four 
sides at the midpoints between the corners, and at the 
center of each supercell, so there is a four-point basis of 
second-layer atoms in each cell. The positions of the first- 
layer atoms can also be seen to be periodically repeated 
in every supercell. 

III. RESULTS 
A. Identification of phases 

Experimentally, there is evidence for liquid, commen- 
surate solid and incommensurate solid phases in the sec- 
ond layer. We now describe the identification of all three 
phases in our simulation. 

To find the liquid phase, we are guided first by the 
torsional oscillator measurements, which detect a liquid 
phase between 0.174 and 0.187 atom/A 2 . We also find 
evidence that densities in this range are liquid in our 
simulation. Figure shows a snapshot of a typical liquid 
density. The second-layer atoms obviously do not possess 
spatial ordering, and the configuration covers the entire 
surface. More direct evidence that the system has a liquid 
phase comes from that static structure factor. Figure ^ 
shows the result of a calculation, which is typical of a 
self-bound liquid, at the coverage 0.1860 atom/A 2 . 

Commensurate and incommensurate solid phases can 
be identified by a similar procedure. A particularly nice 
feature of PIMC is that these solids form on their own, 
without any modifications to the high-temperature den- 
sity matrix, Eq. (^|). In contrast, previous variational 
calculations have used different trial wavefunctions for 
the liquid and solid phases.E3 This can be avoided by 
using a shadow wavefunction, but such calculations have 
not been performed for two-dimensional helium or helium 
films. 

As demonstrated previously,!^ we have observed the 
V / x \fl commensurate solid phase in our simulation 
for temperatures below 1 K. The structure of this phase 
was determined by examining snapshots of the configu- 
rations generated by the simulation. Particle paths of 
the second layer atoms were observed to localize around 
the \fl x \fl lattice sites shown in Fig. [|. We note fur- 
ther that we do not bias the simulation of this solid by 
beginning the configuration at the commensurate lattice 
sites. The existence of the incommensurate solid, which 
occurs at a higher density than the commensurate phase, 
has also been demonstrated. A snapshot of this configu- 
ration generated by our simulation can be found in our 




FIG. 5. Snapshot of a liquid configuration at 0.1778 
atom/A 2 , found using the 5x3 simulation cell with twelve 
active particles and T = 200 mK. Large circles indicate frozen 
first-layer atom sites. The instantaneous configuration of the 
second-layer atoms is indicated by the small circles. 
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FIG. 6. The static structure function for the liquid phase 
at the density 0.1860 atom/A 2 and T = 500 mK with 26 
particles. 
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previous publication. This phase matches the diagram 
shown in Fig. ^. We identify this phase as incommensu- 
rate because no supercell with dimensions less than the 
minimum simulation box dimension can be drawn that 
has both first- and second-layer atoms periodically re- 
peated, in contrast to the commensurate phase. 

The snapshots of the two solids are useful for visualiz- 
ing their structure but are not actual tests for their exis- 
tence. A direct measurement of correlation comes from 
the static structure factor. Results for these calculations 
in the (01) reciprocal lattice direction for the incommen- 
surate and commensurate phases are shown Fig. f?](a) and 
(b). The structure factor is normalized to N act . The lo- 
cations of these peaks give the correct lattice spacings for 
the diagrams shown in the Figs. | and J| The peak for 
the commensurate solid occurs at 1.82 A" 1 , which gives 
the correct lattice constant, 3.99 A, for the y/7 x \fl tri- 
angular solid. Likewise, the peak for the incommensurate 
solid occurs at 1.93 A -1 , corresponding to a lattice con- 
stant of 3.76 A, which is the correct lattice spacing for a 
triangular solid at 0.2083 atom/ A 2 . 



atom/ A 2 , and the layer is uniformly covered by a liquid 
phase from 0.1750 to 0.1905 atom/A 2 . Above this den- 
sity, the liquid coexists with the \fl x v7 commensurate 
solid phase discussed in the previous section. This L-C 
coexistence occurs from 0.1905 to 0.1970 atom/A 2 , and 
is followed by the commensurate phase between 0.1970 
and 0.2032 atom/A 2 . The incommensurate solid phase 
begins to form above 0.2032 atom/A 2 and there is C-IC 
coexistence until 0.2096 atom/A 2 . Above this density, 
until layer promotion to the third layer at 0.212 atom/A 2 , 
the system is completely in the incommensurate phase. 
These results are summarized in Fig. ||(a). 

Before discussing how these ranges were determined, 
we would first like to demonstrate that finite-size effects 
play an unimportant role in the energy values used in the 
Maxwell construction. Figure ||(b) shows the energy per 
particle found using the 8x5 and 5x3 cells. Almost all of 
the points calculated at similar densities in the two cells 
are consistent. The primary "size effect" is the limitation 
on the available densities which may be examined for a 
given simulation cell. 




(a) 
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FIG. 7. The static structure factor calculated in the (01) 
direction for (a) the incommensurate solid at 0.2083 atom/A 2 
and 0.67 K with 32 particles, and (b) the commensurate solid 
at 0.1996 atom/A 2 and 0.50 K with 32 particles. The errors 
are the size of the symbols. 
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FIG. 8. (a) Summary of phase boundaries determined 
from applying the Maxwell construction to the total energy 
of the 8x5 cell. The phases are liquid-gas (L+G), liquid 
(L), liquid-commensurate solid (L+C), commensurate solid 
(C), commensurate solid-incommensurate solid (C+IC), and 
commensurate solid (IC). (b) The energy per particle for the 
5x3 (circles) and 8x5 (squares) cells. 



B. T=0 phase diagram 

Having identified the liquid, commensurate solid, and 
incommensurate solid phases of the second layer, we now 
wish to find the boundaries for each of the phases. We 
are able to identify the following density regions at low 
temperature. At low second-layer coverages, 0.1270 to 
0.1750 atom/A 2 , the system is in a gas-liquid coexistence 
region, which consists of a liquid droplet and a zero den- 
sity gas. The equilibrium density for the liquid is 0.1750 



Phase ranges are determined by using the Maxwell 
double-tangent construction, which identifies unstable re- 
gions associated with the coexistence of two phases. A 
coexistence region at zero temperature in the thermody- 
namic limit will have a total ground state energy that 
is the weighted average of the two constituent phases' 
energy values. In Monte Carlo simulations the energy 
in the coexistence region will lie above the coexistence 
line, either because the system remains in an unphysical 
homogeneous state or because creating the phase bound- 
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ary has a finite energy cost SB We may thus identify a 
coexistence region as the maximum range of densities in 
which all the intermediate energy values lie on or above 
a line connecting the values at the endpoints. We note 
that this version of the Maxwell canstrjiction is somewhat 
different from other applications,Ej'E§E3 which apply the 
Maxwell construction to the free energy dependence on 
atomic area (inverse density) . Our method is appropriate 
for simulations with constant area and varying particle 
number. 

At finite temperatures, the Maxwell construction 
should be applied to the total free energy. Unfortunately, 
the free energy is not directly accessible from the PIMC 
simulation. We instead make use of the fact that at very 
low temperatures the free energy and the energy are ap- 
proximately the same, and become identical at zero tem- 
perature. We can thus apply the Maxwell construction 
to low temperature energy values to determine an effec- 
tively zero temperature phase diagram, provided that the 
values have converged to their zero temperature limits. 
To implement this procedure, we first calculated energy 
values for a range of second-layer densities at 200 mK. 
Selected energy values were recalculated at a higher tem- 
perature, typically 400 mK, and were seen to be within 
error bars of the 200 mK results. This indicates that 
our 200 mK calculations are effectively zero temperature 
results. 
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FIG. 9. The total energy found using the 5x3 simulation 
cell with N act = 8, 9, . . . , 21 and T = 200 mK. The dashed 
line is gas-liquid coexistence line. The solid line indicates a 
coexistence region terminating in an incommensurate solid 
phase. 

The application of the Maxwell construction to the to- 
tal energy values calculated using the 8x5 box has been 
shown in our previous publication.^ Figure |^(a) summa- 
rizes the results. The energy minimum was determined 
to occur at 0.1746 atom/A 2 (30 particles). For compari- 
son, Fig. shows the energy calculations using the 5x3 



cell. These energy values have been shifted by N act e m i n 
for clarity, where e min = -32.754 ± 0.020 K. The en- 
ergy minimum occurs at 0.1778 atom/A 2 (12 particles). 
Note that for both simulation cells the minimum energy 
per particle occurs at nearly the same coverage value, de- 
spite the fact that the 8x5 cell is 2 2/3 times as large 
as the 5x3 cell. In general, we find all the energy values 
calculated with the two cells to be in agreement. See Fig. 

The low density region of the second layer is known 
experimentally to have coexistence between a gas phase 
and a supcrfluid liquid phase. In order to determine the 
gas-liquid coverage range in our simulation, we take the 
gas phase to have zero density at zero temperature and 
thus zero total energy. Two-dimensional calculationsE2l 
confirm that this assumption is correct for low temper- 
atures. A coexistence line can then be drawn between 
the beginning of the second layer, 0.1270 atom/A 2 , and 
the density with the minimum energy per particle, which 
occurs between 0.174 and 0.178 atom/ A 2 in the 8x5 cell. 
The best \ 2 parabolic fit to the energy data around the 
minimum gives pg = 0.1750(6) atom/A 2 for the density of 
minimum energy. The number in parenthesis is the error 
in the last digit. A similar coexistence line can be identi- 
fied in the 5 x 3 cell, Fig. |. We find that p = 0.1752(6), 
so finite-size effects on the liquid density are small. At 
sufficiently low temperatures, this liquid phase will be- 
come superfluid, as will be discussed below. All measured 
energy values for the densities between 0.1270 atom/A 2 
and p lie above the coexistence line, so the system is in 
gas-liquid coexistence for this density range. 

The density of uniform liquid coverage, po, can be com- 
pared to experimental results. ForJT < 0.2 K the second- 
layer heat capacity measurements!!! show a probable gas- 
liquid region roughly between 0.13 and 0.16 atom/A 2 . 
Within the resolution available from the data, this phase 
can terminate anywhere from 0.1600 atom/A 2 up to, but 
not including, 0.1700 atom/A 2 total coverage. Since the 
first-layer coverage in the experiment is between 0.120 
and 0.127 for these densities, gas-liquid coexistence ter- 
minates at second-layer coverages anywhere from 0.033 to 
0.050 atom/A 2 . For comparison, the gas-liquid phase ter- 
minates at the second layer coverage 0.0480(6) atom/A 2 
in our simulation. Superfluidity is first observed in 
the torsional oscillator measurements at 0.174 atom/ A 2 . 
Thus, the superfluid signal, as observed by this technique, 
becomes significant at the coverage where our simulation 
determines that the second layer is uniformly covered by 
the liquid phase. 

The density we determine for uniform liquid coverage 
can also be compared to other simulations. In the two- 
dimensional calculations of Whitlock et al. , the equilib- 
rium liquid coverage is 0.04356 atom/A 2 at zero temper- 
ature. This result is supported tw the low temperature 
results of 2D PIMC calculations.E3 This is slightly below 
our onset coverage, perhaps because we allow for par- 
ticle motion perpendicular to the substrate. Other cal- 
culations for helium films also show a slight increase in 
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the equilibrium density relative to the 2D result. In the 
Monte Casio calculation for the first layer of helium on 
graphite,E3 the equilibrium density is determined to be 
0.0443 atom/A 2 . The effects of wavefunction spreading 
will be even greater in the second helium layer. Wagner 
and Cepecley's simulation of helium adsorbed on solid 
hydrogenc3 also demonstrated that the liquid equilibrium 
density increases when motion perpendicular to the sub- 
strate is allowed. They find a liquid coverage of 0.046(1) 
atom/A 2 , comparable to our result. Thus the calcula- 
tions of films with perpendicular spreading show a trend 
toward higher liquid densities, with the onset density ap- 
proaching the 2D value as the helium-substrate potential 
becomes stronger. From a 2D viewpoint, this can be un- 
derstood as a reduction of the hardcore repulsion, which 
allows for closer crowding. 

At the highest second-layer densities, we can iden- 
tify another unstable region in the total energy values 
of the 8 x 5 cell between 0.2032 and 0.2096 atom/A 2 , 
corresponding to the C-IC mixed phase. As shown 
previouslyjlj the coexistence line can be drawn between 
the total energy values at these two densities. The in- 
termediate energy values lie on or above this line, so the 
region has coexisting phases. In particular, the energy 
value at 0.2080 atom/A 2 was found to lie completely 
above the coexistence line, providing an unambiguous 
signal for coexistence. The range we find is in good 
agreement with the coexistence region 0.2030 to 0.2080 
atom/ A 2 £hat can be determined from the heat capac- 
ity peaksjj This phase coexistence is not a product of 
finite-size effects. The beginning of a similar region may 
be identified between the densities 0.2032 and 0.2117 
atom/A 2 in the 5x3 simulation cell, Fig. ^|. Phase 
coexistence in fact becomes clearer in the 8x5 cell be- 
cause we are able to examine more density values in the 
unstable region. 

The presence of the C phase at 0.1996 atom/A 2 re- 
quires an L-C coexistence region between it and the 
liquid. The region can also be identified in the 8x5 
cell. The endpoints of the L-C phase are 0.1905 and 
0.1969 atom/A 2 . The intermediate energy values lie on 
the coexistence line within error bars. The L-C range 
is in reasonable agreement with the coexistence range 
0.1871 to 0.1970 atom/A 2 determined from heat capap. 
ity measurements. □ Torsional oscillator measurements^!! 
also indicate that the coexistence region begins at about 
0.187 atom/A 2 . The L-C phase cannot be determined in 
the 5x3 cell due to the coarseness of the coverage grid. 



C. Other properties 

Figure |l^ depicts the density profiles for selected layer 
densities. These plots are normalized such that in- 
tegrating p(z) gives the number of particles. Promo- 
tion to the third layer can be clearly observed at the 
highest density shown, 0.2159 atom/A 2 , so we conclude 



that layer promotion occurs between 0.2117 and 0.2159 
atom/ A 2 . This is in excellent agreement with the com- 
pletion density of 0.212 atomiA 2 determined from the 
heat capacity measurements. □ □ A somewhat lower value 
of 0.204 atom/A 2 for third layer promotion is-obtained 
from the isothermal compressibility minimaixEiJ Also of 
note, Whitlock et air 2 ] estimate that promotion to the 
third layer begins at the second-layer coverage of 0.08 
atom/ A 2 , quite close to but somewhat lower than our 
value of 0.085 atom/A 2 . 
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FIG. 10. Density profiles for the second layer found using 
the 5 x 3 cell, with densities 0.1694, 0.1778, 0.1863, 0.1948, 
0.2032, 0.2117, and 0.2159 atom/A 2 . 

The temperature dependence of the energy and su- 
perfluid density at a sample liquid density of 0.1778 
atom/ A 2 have been determined. This coverage corre- 
sponds to a second-layer coverage of 0.0508 atom/A 2 . 
Values were calculated using the 5x3 simulation cell 
with twelve active_particles, and are illustrated in a pre- 
vious publication.E3 The superfluid density is relative to 
the second-layer density. Both the energy and the su- 
perfluid density converge to the ground state for tem- 
peratures below 0.8 K. The slow decay of the superfluid 
density at— higher temperatures is a typical 2D finite- 
size effect .Ell The superfluid density values have been x 2 
fit to the solutioxL, to the Kosterlitz-Thouless (KT) re- 
cursion relationso From the intersection of the fit and 
the KT transition line, we estimate the transition tem- 
perature to be r £c-i~ 0.88K. For comparison, the 2D 
PIMC simulationE] obtains T c = 0.86 ± 0.02if at 0.0508 
atom/ A 2 . 

The specific heat of the liquid, commensurate solid, 
and incommensurate solid phases can be found by dif- 
ferencing the energy per particle with respect to terp-. 
perature. This was shown in our previous publication.t2l 
For the liquid phase, a broad, low peak with a maximum 
value at 1.18 K was found. This is comparable to the ex- 
perimental heat capacity results, □ which have a peak at 1 
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K. For the commensurate solid phase, a specific heat peak 
at about 1.5 K was found. This is comparable to the heat 
capacity measurements^ at similar density values, which 
also have a peak at 1.5 K. This close agreement provides 
some additional evidence that the x V< C phase oc- 
curs in the experiment. Finally, for the IC solid, a peak 
at 0.7 K was obtained, somewhat lower than the peak 
in the heat capacity measurements at the same density, 
which occurs at 1 K. 



IV. SUMMARY 

A number of recent experiments indicate that the sec- 
ond layer of helium on graphite has an interesting phase 
diagram. Torsional oscillator measurements detect-siu. 
perfluidity over a narjcpw density range in this layer. ErEJ 
Neutron scatteringLTu detects an incommensurate solid 
phase at high densities. Heat capacity measurementsO'tl 
have found evidence for liquid-gas coexistence and the 
incommensurate solid phase. The heat capacity data 
also show the existence of an intermediate phase between 
the liquid and incommensurate solid, which is possibly a 
commensurate solid. The existence of this commensu- 
rate solid phase would explain the disappearance of su- 
perfluidity at higher second layer coverages. Motivated 
by these experiments, we have undertaken a simulation 
of this layer. 

In order to study the second layer with Monte Carlo 
for a range of temperatures, it is necessary to develop 
a method that incorporates both particle permutations 
and the effects of the substrate and the solid first layer 
on the second. Permutations are necessary to obtain the 
superfluid phase. The effects on the solid first layer must 
be included since the commensurate second layer solid is 
partially registered with respect to the first layer. First 
layer and substrate effects also play a role in the forma- 
tion of the incommensurate solid phase, which replaces 
the commensurate phase before layer promotion begins. 

We have developed a path integral Monte Carlo 
method that includes the above features. Particle permu- 
tations were included in the. simulation using a method 
developed for bulk heliumEj. We have shown that the 
helium-helium and helium-graphite interactions can be 
incorporated into the simulation by using effective in- 
teractions found from the exact solutions for the inter- 
acting part of the appropriate density matrices. Real- 
istic helium-helium and helium-graphite potentials are 
used to find these effective interactions. For the helium- 
graphite effective interactions, we have shown how this 
solution may be approximated so that off-diagonal ma- 
trix elements may be efficiently and accurately included 
in Monte Carlo sampling. The interaction of the sec- 
ond layer of helium atoms with the solid first layer were 
approximated by placing first layer atoms at triangular 
lattice sites with a lattice spacing that gives the com- 
pleted first layer density. These atoms were located at a 



fixed height above the substrate, given by the minimum 
of the effective helium-graphite interaction. Configura- 
tions of these atoms were not sampled, which allowed us 
to scan second layer densities with a finer grid. There- 
fore, we study the second layer atoms under the influ- 
ence of their mutual interactions and a static potential 
produced by the frozen graphite substrate and the frozen 
first layer helium atoms. This approach ignores effects on 
the second layer from the zero point motion of the first 
layer solid and first layer compression effects. We feel 
this is a reasonable approximation because the relatively, 
high Debye temperature of the completed first layerca 
means that it will be relatively stiff for the temperatures 
of our simulation. Compression effects on the first layer 
by the second are most important for low second layer 
densities^, below the range of our simulation. 

Using the above simulation method, we were able to 
identify, in order of increasing density, superfluid liquid, 
commensurate triangular solid, and incommen- 
surate triangular solid phases from particle configura- 
tions and static structure factor calculations. We also 
calculated the specific heat for each of these phases and 
observed peaks in general agreement with experiment. 

The density ranges at effectively zero temperature of 
each of the second layer phases and their coexistence re- 
gions were determined using the Maxwell construction. 
We found that at low densities, the layer is phase sepa- 
rated into a liquid droplet and a zero density gas. The 
range of this phase is 0.1270 to 0.1750 atom/A 2 . Gas- 
liquid coexistence ends at the equilibrium density for the 
liquid phase. This occurs at 0.1750 atom/A 2 , which is 
the density with the minimum energy per particle. This 
density was found to be insensitive to finite-size effects, 
and is in excellent agreement with the onset of superflu- 
idity determined by torsional oscillator measurements. It 
is also consistent with heat capacity measurements. We 
demonstrated that the liquid phase in our simulation is 
superfluid, and we determined that the transition tem- 
perature was close to the value determined for a purely 
2D superfluid at the same density. 

The helium layer is uniformly covered in our simulation 
by the liquid phase from 0.1750 to 0.1905 atom/A 2 , at 
which point liquid-commensurate solid phase coexistence 
begins. The onset of this coexistence terminates super- 
fluidity, since the growth of the solid phase disrupts the 
connectivity required to detect the superfluid. Exper- 
imentally, liquid-commensurate solid phase coexistence 
has been determined to begin at 0.1870 atom/A 2 by both 
torsional oscillator and heat capacity measurements. We 
determined that the liquid phase is completely replaced 
by the y/7 x y/7 commensurate solid for densities above 
0.1970 atom/A 2 , in good agreement with heat capacity 
measurements. Phase coexistence between the commen- 
surate and incommensurate solid phases begins at 0.2032 
atom/ A 2 . For coverages above 0.2080 atom/ A 2 , the in- 
commensurate solid is the only phase occurring until 
layer promotion. These ranges for the solid coexistence 
and the incommensurate solid are in agreement with the 
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heat capacity measurements. The density ranges for all 
the second layer phases described above are summarized 
in Fig. ||(a). Finally, we observed layer promotion for 
coverages above 0.2117 atom/ A 2 , in excellent agreement 
with experiment. 
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